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Abstract. We consider a stochastic environment with two time scales and outline a general 
theory that compares two methods to reduce the dimension of the original system. The first method 
involves the computation of the underlying deterministic center manifold followed by a "nai've" 
replacement of the stochastic term. The second method allows one to more accurately describe the 
stochastic effects and involves the derivation of a normal form coordinate transform that is used 
to find the stochastic center manifold. The results of both methods are used along with the path 
integral formalism of large fluctuation theory to predict the escape rate from one basin of attraction 
to another. The general theory is applied to the example of a surface flow described by a generic, 
singularly perturbed, damped, nonlinear oscillator with additive, Gaussian noise. We show how both 
nonlinear reduction methods compare in escape rate scaling. Additionally, the center manifolds are 
shown to predict high prc-history probability regions of escape. The theoretical results are confirmed 
using numerical computation of the mean escape time and escape prehistory, and we briefly discuss 
the extension of the theory to stochastic control. 
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1. Introduction. It has long been known that noise can have a significant effect 
on deterministic dynamical systems. For example, given an initial state starting in 
some basin of attraction (defined as the set of initial conditions from which the system 
approaches a corresponding locally stable attractor as time evolves to infinity), noise 
can cause the initial state to cross the basin boundary and move into another, distinct 
basin of attraction p i [14 ] |39 ] [36 j [35] . 

There are several points of view one might consider when investigating the effect 
of noise on a dynamical system, including stochastic resonance [24] and finite noise 
effects [5]. In this article, we consider yet another point of view, namely the effect of 
arbitrarily small noise on the escape of a particle from a potential well. In this case, 
one can apply large fluctuation theory [20] [12l [HI [36] . 

Many of the underlying deterministic systems found in [12] [HI [39] [361 ES] have 
parameter regimes in which multiple attractors give rise to noise-induced escape from 
one attractor to another. Such systems may be analyzed globally by considering the 
Hamiltonian theory of large fluctuations or by considering escape from attracting 
potential wells along most probable exit paths [23] [26] [37] [28] . 

Through the use of a path integral coupled with variational methods, it is possible 
to compute the probability densities of the trajectories of the system. In particular, 
for sufficiently small noise, one can find the trajectories which escape from a basin 
of attraction due to stochastic effects. The most probable escape trajectory is the 
optimal escape path of a state residing in a basin of attraction. 

Many researchers have investigated how noise affects physical and biological phe- 
nomena, including lasers [8] [32] [27] , epidemics and control [4] [17l [331 US] j ana - neu- 
rons [30] , Yet another important application in many fields is that of sensing in 
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stochastic environments. Improved environmental sensing and prediction can be 
achieved through the incorporation of continuous monitoring of the region of interest. 
For example, one could monitor the stochastic ocean using autonomous underwater 
gliders [48l [47[ [19] . However, to do this, one must understand both the dynamics and 
control of the gliders. 

Extending the lifetime (energy optimization problem) of sensing devices (e.g. 
gliders) in stochastic environments such as the ocean requires an understanding of the 
effect of the environmental forces on both the devices and the region being monitored. 
The ocean dynamics are high-dimensional and stochastic. Therefore, as a first step 
towards using the underlying ocean structure to optimize a sensor's energy usage, 
we will outline a general theory that provides two methods to obtain a reduction in 
the dimension of the stochastic system. The manifold equations that are found using 
these methods can then be used to determine the optimal escape path and escape 
rate. 

Our formulation uses large fluctuation theory [20] to determine the first passage 
times in a multi-scale environment. For a vector field that has relaxation times on the 
same scale, it is clear how to use the theory to generate an optimal path of escape, and 
this theory has been applied to a variety of Hamiltonian and Lagrangian variational 
problems [3D [H] - 

However, technical issues may arise when one wishes to determine the projection 
of noise that is needed to perturb the dynamics restricted to the lower-dimensional 
manifold. To address this, several approaches have been developed to understand 
dimension reduction in systems that have well separated time scales. For a system with 
certain spectral requirements, the existence of a stochastic center manifold was proven 
in [B]. Non-rigorous stochastic normal form analysis (which leads to the stochastic 
center manifold) was performed in [34] [9[ [40l [41] . Rigorous theoretical analysis of 
normal form coordinate transformations for stochastic center manifold reduction was 
developed in [21[T]. Later, an alternative method of stochastic normal form reduction 
was developed [42] , in which any anticipatory convolutions (integrals into the future of 
the noise processes) that appeared in the slow modes were removed. Since this latter 
analysis makes the construction of the stochastic normal form coordinate transform 
more transparent, we use this method to derive the reduced stochastic center manifold 
equation. 

The layout of the paper is as follows. The general theory of deterministic and 
stochastic center manifold reduction is described in Sec. [51 The first method used 
to reduce the dimension of the system involves the derivation of the center manifold 
equation [7] of the associated deterministic system followed by the "nai've" replacement 
of the stochastic term. The second method, which allows one to more accurately 
describe the effect of the noise, involves the derivation of a normal form coordinate 
transform [55] that is used to find the stochastic center manifold equation. Section [5] 
also describes how the two center manifolds resulting from the two methods can be 
used along with the theory of large fluctuations to analytically find the optimal escape 
path of the particle along with its escape rate. The general theory of Sec.[2jis applied 
to a specific example given by a singularly perturbed, damped, Duffing oscillator with 
additive, Gaussian noise in Sec. [3[ This section contains analytical results and their 
comparison with numerical computation. The conclusions are contained in Sec.UJ 

2. General Theory. We consider the following general (m + Tridimensional 
system of stochastic differential equations with two well-separated time scales: 
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(2.1a) x = Ax + F(x,y,$), 

(2.1b) ey = By + G(x,y,*), 

where e is a small parameter, x(i) g R m , y(t) £ M. n , $(t) and *&(t) describe stochastic 
forces with adjustable intensity, A and B are constant matrices, and F and G are 
stochastic, nonlinear functions. 

2.1. Deterministic Center Manifold. To begin, we remove the stochastic 
terms from Eqs. (|2~Ta| and (|2.1b|) so that F = F(x, y) and G = G(x, y). Let t = er. 
Denoting ' as d/dt and ' as d/dr, then the deterministic form of Eqs. (|2.1a|) and (I2.1bl) 
is transformed to the following system of equations: 

(2.2a) x' = e [Ax + F(x, y)] , 

(2.2b) y' = By + G(x,y), 

(2.2c) e = 0. 

To recast the problem in a more general framework, we treat e as a state variable, 
let A = eA and F = eF, and write Eqs. (|2~2aj) - ([2~2c)) as 

(2.3a) x' = Ax + F(x,y,e), 

(2.3b) y' = By + G(x,y), 

(2.3c) e = 0. 

If A and B arc constant matrices such that all of the eigenvalues of A have zero real 
parts, while all of the eigenvalues of B have negative real parts, then the system will 
rapidly collapse onto a lower-dimensional manifold given by center manifold theory [7] . 
Furthermore, we will consider examples where the solution decays throughout the 
transient and then stays close to the lower-dimensional manifold. 
If the center manifold is given by 

(2.4) y = h(x,e), 

then substitution of Eq. (|2.4[) into Eq. (|2.3b[) leads to the following center manifold 
condition: 

(2.5) h x [Ax + F(x, h(x, e), e)] = Bh(x, e) + G(x, h(x, e)), 

where h x denotes the partial derivative of h with respect to x. Although it is generally 
not possible to solve Eq. (|2 . 5[) for h, one can approximate the center manifold by 
expanding h in the following way: 

(2.6) h(x, e) = ho(x) + ehx(x) + e 2 h 2 (x) + 0(e 3 ). 

Typically, this approximation of h(x, e) is found by substituting Eq. (|2.6p into the 
center manifold condition [Eq. (|2.5p ] and matching coefficients. 
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2.2. Optimal Escape Path and Escape Rate. Starting with the center man- 
ifold equation given by Eq. (|2.6[) to a particular order, the dynamics on the center 
manifold can be found by substitution of Eq. (|2.6[) into Eq. (|2.2ap [using the relation 
given by Eq. (|2.4[) ]. Therefore, the dynamics are determined by 

(2.7) x' = e [Ax + F(x, h(x, e))] = eH(x, e), 
and use of the relation between t and r leads to the following: 

(2.8) x = Ax + F(x,h(x,e)) = H(x,e). 

Equation (|2.8|) is a deterministic equation. However, now that we have reduced 
the dimension of the problem, we return to considering a stochastic problem by 
"naively" adding a noise vector to the right-hand side of Eq. (|2.8|) so that one has 

(2.9) x = H(x,e) + \/2£>#(t), 

where D is the noise intensity. Each of the noise components, cf>i, of $ in Eq. (|2.9I) 
describes a stochastic white force that is characterized by the following correlation 
functions: 

(2.10a) {4>i(t)) = 0, 



(2.10b) (&(*)&(*0) =S(t-t')S ij . 

The following analysis to determine the optimal escape path may be performed 
using Eqs. (|2.9p - (|2.10b|) , a system with an arbitrary number of degrees of freedom [15] . 
However, since we ultimately are interested in applying this general theory to a 
two-dimensional (2D) surface flow that reduces via the center manifold to a one- 
dimensional (ID) equation, for simplicity we consider the ID version of Eq. (|2.9D 
given as 

(2.11) x = H(x,e)+V2D4>(t), 

where (j)(t) is characterized by the correlation functions given by Eqs. (|2.10ap and 

(l2~Tobl) . 

We assume that H(x,e) is associated with a potential function U(x, e), 
[H(x,e) = —dU(x,e)/dx] with stable states (attractors) located at x = x a and an 
unstable state (saddle) located at x = x s . Then Eq. (|2.1ip corresponds to a Langcvin 
equation of a particle in an over-damped potential well. 

In the absence of noise, a particle located in the potential well will approach the 
stable, attracting state. However, the noise may organize as an effective force which 
acts on the particle and "pushes" the particle from the attractor to the saddle located 
at the top of the potential well barrier. The path along which the particle leaves the 
basin of attraction due to such an effective noise force is an escape path. 

Given that 4>(t) is uncorrelated Gaussian noise, the probability of optimal escape 

is 



(2.12) P[x csc ] = Cexp / <^ pt dt \ = Ccxp (-R/D) 
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(2.13) R=- J 4> 2 opt dt = ^ J L(x,x;t)dt, 



and 4> op t denotes the stochastic fluctuations corresponding to the trajectory that 
moves along the optimal escape path. In Eq. (|2.12|) . C is a pre-factor that depends on 
the noise intensity. In order to maximize the probability of escape, one must minimize 
the exponent given by Eq. (|2.13[) . Using the Euler-Lagrange equation of motion, it is 
now possible to solve for the optimal escape path. 

It also is possible to derive an expression for the escape rate from an attractor to 
the saddle, which is located at the top of the potential well barrier. In general, if the 
stochastic trajectory of a particle is given by 

(2.14) X = - d -^ + V2Dm, 

with U (x, e) some potential (as we have assumed), then for Gaussian noise, the escape 
rate from the attractor (located at x = x a ) to the unstable saddle (located at x = x s 
at the top of the barrier) is given by 



(2.15) W(D) = V\Q(xs,e)\Q(x a i e)_ ^ r AU/L) . 

where Q(x,e) = d 2 U(x,e)/dx 2 and AU is the activation energy of escape (depth of 
the potential well) [53] . 

If W{D) is the escape rate, then 1/W(D) gives the mean escape time, and the 
natural log of the mean escape time is therefore 

2?r \ AU 



(2.16) log — — = log 



W(D)J *\^\Q(x s ,e)\Q(x a ,e)J D' 

It should be noted that to ensure the particle escapes from the potential well, 
one could compute the escape rate (and the mean escape time) from the attractor 
to a point located somewhere in the second potential well (i.e. the particle escapes 
from the first basin of attraction if it climbs out of the potential well to the saddle 
located at the top of the potential well barrier and then continues past the saddle to 
some point located in the second basin of attraction). However, this movement of the 
point at which one claims the particle has escaped from the potential well will create 
a change in the pre-factor of Eq. (|2.15[) along with a corresponding change in the first 
term on the right-hand side of Eq. (|2.16[) [25] . 

As we have previously noted, we consider systems whose solution decays exponen- 
tially throughout the transient and then stays close to the lower-dimensional center 
manifold. There are no secular terms in the asymptotic expansion since we are not 
looking at periodic orbits, and the result is valid for all time. Moreover, any noise drift 
on the center manifold will result in bounded solutions due to sufficient dissipation 
transverse to the manifold. This behavior is in direct contrast to the finite-time solu- 
tions on manifolds of relaxation-type oscillators whereby the time scale of escape must 
be shorter than the lifetime of the trajectories on the manifold [I0|[22]. This lifetime 
issue does not pertain to the systems we consider since there will be no oscillations. 
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2.3. Stochastic Center Manifold and the Normal Form Coordinate 
Transform. In a manner similar to that shown in Sec. 12.11 the stochastic system 
given by Eqs. (|2.1ap and (|2.1b[) can be transformed [3] to the form given by Eqs. (|2.3a|l - 
(|2.3cj) . where now F = F(x, y, e, 3>) and G = G(x, y, If A and B satisfy the same 
spectral conditions as for the deterministic system, and if the stochastic time depen- 
dence found in F and G is due to independent white noise processes, then there exists 
a stochastic center manifold for the original stochastic system [6]. 

One method for computing the stochastic center manifold for systems with both 
fast and slow dynamics uses the construction of a normal form coordinate transform 
that not only reduces the dimension of the dynamics, but also separates all of the fast 
processes from all of the slow processes [42]. While this type of normal form coordi- 
nate transform may be used to find deterministic center manifolds, the application of 
this transform to stochastic systems is particularly interesting since white noise has 
fluctuations on all scales. 

There are many publications, such as pHl SOI E] which deal with the simplifica- 
tion of a stochastic dynamical system using a stochastic normal form transformation. 
In these articles, the noise term is multiplied by a small parameter, and therefore, 
the resulting stochastic normal form is a perturbation of the deterministic normal 
form. Furthermore, one can find in 141] normal form transformations that involve 
anticipative noise processes. However, these integrals of the noise process into the 
future were not dealt with rigorously in [9l [41] . 

Rigorous, theoretical analysis to support normal form coordinate transforms (and 
center manifold reduction) was developed in [21 [1] . In this work, the technical prob- 
lem of the anticipative noise integrals also was dealt with rigorously. Later, another 
stochastic normal form transformation was developed [42] - This new method is such 
that "anticipation can ... always [be] removed from the slow modes with the result 
that no anticipation is required after the fast transients decay" (Ref. [42], pp. 13). An 
advantage of removing anticipation is the simplification of the normal form. Nonethe- 
less, this simpler normal form retains its accuracy with the original stochastic system. 
Furthermore, when modeling the macroscopic behavior of microscopic, stochastic sys- 
tems, it is desirable to avoid anticipation in the normal form |42j . 

It is important to note that the normal form is valid for all time since it is just a 
coordinate transform. Furthermore, the dynamics also are valid for all time as long 
as the truncation error is small enough for the problem of interest. 

In the example of Sec. El we shall use the method of [42] to simplify our stochas- 
tic dynamical system to one that emulates the long-term dynamics of the original, 
multiplc-timc-scale system. The method involves five principles, which we recapitulate 
here for the purpose of clarity. The principles are as follows: 

1. Avoid unbounded, secular terms in both the transformation and the evolution 
equations to ensure a uniform asymptotic approximation. 

2. Decouple all of the slow processes from the fast processes to ensure a valid 
long-term model. 

3. Insist that the stochastic slow manifold is precisely the transformed fast pro- 
cesses coordinate being equal to zero. 

4. To simplify matters, eliminate as many as possible of the terms in the evo- 
lution equations. 

5. Try to remove all fast processes from the slow processes by avoiding as much 
as possible the fast time memory integrals in the evolution equations. 

In practice, the original stochastic system of equations (which satisfy the necessary 
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spectral requirements) in (x, y) coordinates is transformed to a new (X, Y) coordinate 
system using a stochastic coordinate transform as follows: 

(2.17a) x = X + £(X,Y,t), 

(2.17b) y = Y + T7(X,Y,i), 

where the specific form of Eqs. (|2. 17a|) and (|2.17b[) is chosen to simplify the origi- 
nal system according to the five principles listed previously. The terms £(X,Y, t) 
and r/(X, Y, t) are found using an iterative procedure that will be demonstrated us- 
ing the singularly perturbed, damped, stochastic Duffing oscillator model in Sec. [3] 
Theoretical details can be found in [12] . 

3. Example - Singularly Perturbed Stochastic Duffing Oscillator. We 

consider the following singularly perturbed, damped, Duffing oscillator system with 
additive noise: 

(3.1a) x = y + y/2D<t>{t), 

(3.1b) ey = (x — x 3 — y), 

where D is the noise intensity and <f>(t) describes a stochastic white force that is 
characterized by the correlation functions given in Eqs. (|2.10a[) and (|2.10bl) . 

In this example, the noise is added only to the x equation. Additionally, one 
could consider two other scenarios. In the first scenario, noise is added to both the x 
and y equations, while in the second scenario, noise is added only to the y equation. 
To implement the first scenario, one adds y/2D(f>i(t) to the right-hand side of the x 
equation and V2D(f>2 (t) to the right-hand side of the y equation [where (j>i(t) and <f>2(t) 
describe stochastic white forces of intensity D that are characterized by the correlation 
functions given in Eqs. (|2.10a|) and (|2.10b|) ]. To implement the second scenario, noise 
is added only to the y equation. However, unlike the example [Eqs. ([3. La[) and ([3. lb)) ] 
and the first scenario, in this case the noise term must be scaled by y/e and the 
potential function must be scaled by e [44] . Although the following results pertain to 
Eqs. (|3.1aj) and (|3.1b[) . we have checked that one obtains similar results for the other 
two scenarios. 

The system given by Eqs. (|3~Ta| and ([3~Tb]) is very strongly damped when eCl. 
For the case of an under-damped system, one should consider the dynamics in the 
limit of weak damping [17] . 

3.1. Deterministic Center Manifold. Following the general theory of Sec. 
12.11 we consider the deterministic form of Eqs. (|3.1a[) and (|3.1b[) by setting (f)(t) = 0. 
The slow manifold is found by setting e = in Eq. (|3.1b[) . Solving for y gives the 
equation of the slow manifold as y — x — x 3 [which corresponds to Hq(x) in Eq. (|2.6p ]. 
Substitution of this into the deterministic form of Eq. (|3.1a[) gives the dynamics along 
the slow manifold as x = x — x 3 . 

If, as in Sec. 12. 11 we let t — er and denote ' as d/dt and ' as d/dr, then Eqs. (|3.1ap 
and (|3.1b[) (with cf)(t) = 0) are transformed to the following system: 

(3.2a) x 1 = ey, 



(3.2b) 



y' = x- x 3 - y, 
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(3.2c) e' = 0. 

Rearrangement of Eqs. (|3.2a|) - (|3.2c|) leads to a system described by constant matrices 
A and B that satisfy the spectral requirements of Sec. 12.11 Furthermore, since the x 
and e variables are associated with the A matrix (eigenvalues with zero real parts), 
and the y variable is associated with the B matrix (eigenvalues with negative real 
parts), we know that the center manifold is given by y = h(x, e). 

The center manifold condition is given by Eq. (|2.5[) . and we approximate the 
center manifold [Eq. (|2 . 6[) ] as follows: 

(3.3a) h{x, e) = h (x) + ehi(x) + e 2 h 2 (x) + 0{e 3 ) 



= c + c ie + c w x + c 02 e 2 + c n xe + c 20 x 2 
(3.3b) + c 03 e 3 + c 12 xe 2 + c 21 x 2 e + c 30 a; 3 + 0(j 4 ), 

where Co, Coi, cio, cq 2 , ■ ■ ■ are unknown coefficients, and 7 = \(x, e)| so that 7 provides 
a count of the number of x and e factors in any one term. The center manifold 
condition for this example is given by 

(3.4) £) [eh(x, e)] = -h{x, e) + x-x 3 . 

By substituting Eq. (|3.3b[) into Eq. (|3.4p and matching the different orders to 
find the coefficients, one finds the following center manifold equation (expanded to 
sixth-order): 

h(x, e) = x - ex + 2e 2 x - x 3 ~ 5e 3 x + Aex 3 + 14e 4 x 
(3.5a) - 20e 2 x 3 - 42e 5 x + I04e 3 x 3 - 3ex 5 + C( 7 7 ) 

= x - x 3 + e(-x + 4x 3 - 3a; 5 ) + e 2 (2x - 20a; 3 ) 
(3.5b) + e 3 (-5a; + 104a; 3 ) + e 4 (14a;) + e 5 (-42a;) + C( 7 7 ). 

Note that by letting e = 0, one recovers the zero-order approximation, ho(x) (the slow 
manifold). In addition, since e is now a state variable, the first nontrivial correction 
term to the zero-order approximation is a quadratic term. 

3.2. Optimal Escape Path and Escape Rate. Consider the third-order cen- 
ter manifold equation given by h(x, e) = x — x 3 — ex + 2e 2 x [see Eq. (|3.5a[) ]. Following 
Scc. 12.21 the dynamics on the center manifold are given by x' — eh(x, e) [see Eq. (|2.7[) ]. 
Use of the relation t = er leads to the following deterministic equation: 

(3.6) x = h(x, e) = x — x 3 — ex + 2e 2 x. 

We "naively" add the noise term cj)(t) to the right-hand side of Eq. (|3.6|) as in 
Eqs. ([279]) and (f2~TTj) so that 

(3.7) ± = a;-a; 3 -ex + 2e 2 3; + v / 2Dc/'(t). 

Equation (|3.7|) corresponds to a Langevin equation of a particle in an over- 
damped quartic potential well with stable states (attractors) located at x = x a = 
±\/l — e + 2e 2 and an unstable state (saddle) located at x = x s = 0. 
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The probability of optimal escape is given by Eqs. (|2.12[) and (|2.13[) . Solution of 
the Euler-Lagrangc equation of motion leads to the following optimal escape path: 



(3.8) 



.4 



1 + 3 exp [2 At) ' 



where A = 1 — e + 2e 2 . Note that a; esc — > \J~A as t — > — oo while x csc — ► as t — ► oo. 
This path is a heteroclinic orbit from the basin of attraction located at x = \J~A to 
the saddle located at x = 0. 
Using Eq. (|3.8|) along with 

(3.9) ^csc — -^csc ^csc e ^csc ~t~ 2c Xcsc ~t~ V^2Z?0 O pt (t) , 

we find that the optimal noise is given by the following: 

-6Aexp(2At) 



(3.10) opt (t) 




3 exp {2 At) 



1 + 3 exp (2 At) 



Since A is a function of e, the shape of </> pt(*) w iU be affected by the value of e. With 
D = 0.05, Fig. I5.1f a) shows <j> op t(t) for various values of e. To obtain a clearer view, 
Fig. I5.1f b) shows a section of Fig. IQT a). One can see in Figs. IS-ll fa) and !5.1f b) that 
e has an effect on both the pulse width and amplitude. Starting with e = 0.02, the 
pulse amplitude decreases monotonically and the pulse width increases monotonically 
as e increases to e = 0.25 (the value of e for which A is minimized). As e continues to 
increase beyond e = 0.25, the pulse amplitude increases monotonically and the pulse 
width decreases monotonically. 

Using the theory outlined in Sec. 12. 21 we now derive an expression for the escape 
rate from one of the attractors to the saddle in order to predict the change in escape 
rate caused by varying the value of e and D. 

Since the stochastic ( "naive" ) , third-order center manifold dynamical equation 
given by Eq. (|3 . 7|) has the form of Eq. (|2.14[) . then the escape rate from the attractor 
located at x — x a = yl — e + 2e 2 to the saddle located at x = x s = can be found 
using Eq. (|2.15[) . The escape rate is given as follows: 



(3.11a) W(e, D) = VI - 1 + ^ - 2^|(2 - 2e+^) y ^ 

2tt 



(3.11b) 



AU = 



-l + 2e-5e 2 + 4e 3 -4e 4 



Appendix IA1 contains similar expressions for the escape rate that are found using the 
fourth-order and fifth-order stochastic ("naive") center manifold dynamical equations. 

3.3. Stochastic Center Manifold and the Normal Form Coordinate 
Transform. To more accurately describe the stochastic effects, we will derive the 
normal form coordinate transform (and thus the stochastic center manifold) for the 
singularly perturbed, stochastic Duffing system given by Eqs. (|3.1a[) and (|3.1b[) . As 
demonstrated previously, use of the t = er transformation leads to the following 
system: 
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(3.12a) x' = e(y+ V2D(j>) = e(y + cr0), 

(3.12b) y' = x~x 3 -y, 

(3.12c) e' = 0, 

where a is the standard deviation of the noise intensity D = a 2 /2. 

The construction of the normal form is quite tedious and complicated. However, 
the result allows one to determine if there are any noise terms that cause a signifi- 
cant difference between the average stochastic center manifold (the stochastic center 
manifold generally fluctuates about an average location) and the deterministic center 
manifold. 

For this problem, it turns out that the noise terms that could lead to a difference 
between the deterministic and average stochastic center manifolds occur at very high 
order in the normal form expansion. Therefore, the correction to the deterministic 
center manifold is minimal, and we expect that the deterministic results of Sec. 13. H and 
Sec. 13.21 will agree very well with numerical computations using the original stochastic 
system [Eqs. (|3Ta|) and (|3.1bj) ]. 

We proceed by showing how to use the method of [42] described in Sec. 12.31 
to construct a normal form coordinate transform that separates the slow and fast 
dynamics of Eqs. (|3.12a[) and (|3.12b[) . In what follows, we outline the steps involved 
in the first iteration, while details regarding the higher iterations are provided in the 
appendices. 

3.3.1. First Iteration. We begin by letting 
(3.13a) x ps X, 

(3.13b) X' ps 0, 

and by finding a change to the y coordinate (fast process) with the form 
(3.14a) y = Y + r)(T,X,Y) + ..., 



(3.14b) Y' = -Y + G(t,X,Y) + ..., 

where 77 and G are small corrections to the coordinate transform and the corresponding 
evolution equation. Substitution of Eqs. (|3.13a[) - (|3.14b|) into Eq. (|3.12b|) gives the 
following equation: 

<"»> ^IMf+l^-- 3 

Replac ing Y' = dY/dr with —Y + G [Eq. ifsTlibJ) ]. noting that dX/dr = [Eq. 
(I3.13b[) ]. and ignoring the term drj/dY ■ G since it is a product of small corrections 
leads to the following: 



(3.16) 
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Equation (|3.16|) must now be solved for G and r\. In order to keep the evolution 
equation [Eq. (|3.14b[) ] as simple as possible (principle (0| of Sec. 12. 3p . we let G = 0, 
which means that the coordinate transform [Eq. (|3.14a[) ] is modified by 77 = X — X 3 . 
Therefore, the new approximation of the coordinate transform and its dynamics are 
given by 

(3.17a) y = Y + X - X 3 + 0{( 2 ), 

(3.17b) Y' = -Y + C(C 2 ), 

where £ = \(X, Y, e, cr)| so that ( provides a count of the number of X, Y, e, and a 
factors in any one term. 

3.3.2. Higher Iterations. The construction of the normal form continues by 
seeking corrections, £ and F, to the x coordinate transform and the X evolution using 
the updated residual of the x equation [Eq. (|3.12ap ]. and by seeking corrections, rj 
and G, to the y coordinate transform and the Y evolution equation using the updated 
residual of the y equation [Eq. (|3.12bj) ]. Details regarding the second iteration can be 
found in Appendix IbI 

The derivation of £ and F in the second and fourth iterations along with the 
derivation of 77 and G in the third iteration leads to the following updated approx- 
imation of the coordinate transforms and their corresponding evolution equations: 

y =Y + X - X 3 + e (-X + AX 3 - 3X 5 ) 

+ ecr (-e- T * cf) + 3X 2 e~ T * 0) + 3e 2 XY 2 + C(C 3 ), 

Y' = - Y + e (-Y + 3X 2 Y) + C(C 3 ), 



(3.18a) 
(3.18b) 



x =X - eY + ( 2 (Y - 'iX 2 Y) 
(3.19a) + e 2 a (e~ T * cf> - 3X V T * (f) + C(C 4 ), 

X' =e (X - X 3 ) + ea<j) + e 2 (-X + AX 3 - 3X 5 ) 
(3.19b) + e 2 a (-0 + 3X 2 <jy) + G(C 4 ), 

where 

T 

(3.20) e~ T *(/>= J exp[-(T-s)}ct>(s)ds, 

—00 

Details regarding the derivation of Eqs. (|3.18ap and (|3.18b[) can be found in Ap- 
pendix [Cl while details of the derivation of Eqs. (|3.19ap and ()3. 19b[) can be found in 
Appendix [D] 

One can continue this iterative procedure to obtain higher order terms in the 
expansions of the coordinate transform and normal form. For the stochastic Duffing 
system under consideration, the fifth and sixth iterations lead to updated approxi- 
mations of the x and y coordinate transforms (along with their associated evolution 
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equations) that are extremely long and complicated. These approximations can be 
found in Appendix El 

In the higher order transform given by Eqs. (|E.la[) - (|E.ld|) . one can see the 
appearance of quadratic noise terms. For example, one can see terms of the form 
e~ r * {eT T * <j>) 2 in the coordinate transforms [Eqs. (jE.laJ) and (|E.lc[) ]. and one can 
see terms of the form cf>e~ T * <f> in one of the evolution equations [Eq. (|E.ld[) ]. This 
quadratic noise is important because it leads to the creation of a deterministic drift 
within the slow dynamics [321 HI] • Furthermore, the stochastic center manifold gener- 
ally undergoes fluctuations about a mean or average location. This average stochastic 
center manifold is usually different from the deterministic center manifold, and it is 
the quadratic noise process that generates this difference. 

3.3.3. Comparison with Deterministic Center Manifold and Effect of 
Quadratic Noise. Letting Y = and a — in Eqs. (|E.la|) and (|E.lc[) leads to the 
following deterministic center manifold equation: 

(3.21a) x =X, 

y =X - X 3 + e(-X + 4a; 3 - 3X 5 ) + e 2 (2X - 20X 3 + 42X 5 ) 
(3.21b) + e 3 (-X + 16X 3 - 66X 5 + 96X 7 - 45X 9 ) + 0{e 3 ). 

Comparison of Eqs. (|3.21ap and (|3.21b[) with Eq. (|3.5b[) shows agreement through the 
0(e 2 ) terms. There appears to be a discrepancy at order 0(e 3 ). However, we have 
checked that this apparent discrepancy is resolved by expanding the stochastic normal 
form coordinate transform to even higher order. For example, the seventh iteration 
will yield a — Ae 3 X term in the y coordinate transform. When added to the existing 
— e 3 X term, there is an agreement with the — 5e 3 a; term in Eq. (|3.5b|) . 

Letting Y = in Eqs. (jE.laj) and (jE.lcj) leads to the stochastic center manifold 
equation. If one takes the expectation of this stochastic center manifold equation and 
uses the following identities [42] : 

(3.22) E[e ±T *0] =e ±r * E[<j>], 

(3.23) J B[(e ± ^0) 2 ] = i, 
then one obtains the following: 

E[y] =X - X 3 + e(-X + AX 3 - 3X 5 ) + e 2 (2X - 20X 3 + 42X 5 ) 
+ e 3 (-X + 16X 3 - 66X 5 + 96X 7 - 45X 9 ) 

(3.24) +e i a 2 (-3X/2 + 9X 3 -27X 5 /2), 

where the 0(e 4 a 2 ) terms are associated with the quadratic noise terms in Eq. (jE.lap . 

The average stochastic center manifold equation given by Eq. (|3.24[) can now be 
used in conjunction with the theory outlined in Sec. l2.2l to find an analytical expression 
for the escape rate. This calculation has been performed, but because the noise effects 
occur at such high order [0(e 4 a 2 )], the correction to the stochastic ("naive") result 
is minimal. 
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3.4. Numerical Computation of Escape Time. When a = (no noise), 
the original, singularly perturbed problem given by Eqs. (|3.1ap and (|3.1b[) . has three 
equilibrium points given by (—1,0), (0,0), and (1,0). At the initial time, t = 0, a 
particle is randomly placed near the stable, attracting point (1,0) within a circle of 
radius 0.1 centered at (1,0). 

Equations (13. lap and (|3.1b[) are numerically integrated using a stochastic fourth- 
order Runge-Kutta scheme [531 ES] with a constant time step size, St, that depends on 
the value of e (St = 0.01 for e < 0.1, while St = 0.1 for e > 0.1), and the time needed 
for the particle to escape from the basin of attraction is determined. This escape time 
is based on either the time it takes the particle to cross the x < —0.2 barrier, which 
means the particle has escaped across the unstable saddle, and has entered the second 
basin of attraction with stable, attracting point (—1, 0), or when the maximum time 
(10, 000, 000 for St = 0.01 and 100, 000, 000 for St = 0.1) has been reached. 

This computation was performed for 10, 000 particles, and the mean escape time 
was determined. In these computations, e ranged from e = 0.02 to e — 1.5, and a 
ranged from a = 0.26 to a = 0.5. Figure [S~2l shows a contour plot of the natural log 
of the mean escape time plotted for the above range of e and 1/D = 2/cr 2 values. 

By taking a vertical slice of Fig. 15. 2[ one can look at a plot of the natural log of 
the mean escape time versus 1/D for a fixed value of e. Figure |5~3T a) shows a vertical 
slice of Fig. l5.2l taken at e = 0.1 along with a line of best fit through these numerically 
computed data points. The slope of the best fit line is m = 0.2505. From Eq. (|2.16p . 
we see that if one plots the natural log of the mean escape time vs. 1/D, then the 
theoretical slope of this line is given by AU, the depth of the potential well. If e = 0, 
the depth of the potential well corresponding to Eqs. (|3.1ap and (|3.1bp is AU = 0.25, 
which compares very well with the numerical result for e = 0.1. 

Figure [5~3T b) shows several vertical slices of Fig. 15.21 taken at various values of e. 
As in Fig. I5.3f a). there is a line of best fit through the data points corresponding to 
each choice of e. 

Numerical computations have been performed for values of e as small as 0.02. 
The slopes of the lines of best fit through the data for these small values of e are very 
close to 0.25. The data and their lines of best fit are not shown in Fig. I5.3f b) since 
the plots would obscure one another (and they would all be obscured by the e = 0.1 
plot). 

One can see in Fig. I5.3f b) that as the value of e increases (which means that the 
system moves further and further away from the over-damped regime) , the slope of the 
line of best fit decreases. The slope values are as follows: for e = 0.2, m = 0.2466; for 
e = 0.3, m = 0.2398; for e = 0.4, m = 0.2287; for e = 0.5, m = 0.2165; and for e = 1.0, 
m = 0.1549. When the system is in an under-damped regime, the slope of the best 
fit line no longer agrees with the theoretical value when e = (m = 0.25). However, 
there still is a nice scaling behavior. One should note that in general the escape rate 
computed using the one-dimensional system found by letting e = may be different 
than the escape rate of the associated two-dimensional system with e — > [11]. A 
specific example may be found in the case of extinction processes |46j , where the rate 
to extinction has no limit when a singular parameter becomes small. 

As shown in Sec. [2^1 Eqs. (|3.11a[) and (j3.11bj) can be used to find analytical 
values of the natural log of the mean escape time. By varying the value of D and 
by fixing the value of e, we can compare the analytical mean escape time found 
using Eqs. (|3.11ap and ()3.11b[) [or Eqs. (|A~2a)) - (IA.2dP ; Eqs. (|A~4a) - (|A.4dp ] with the 
numerically computed mean escape time shown in Figs. [3~2l and [5~3l 
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With e = 0.02, Fig. 15.4( a) shows a comparison between the numerically computed 
mean escape time (Fig. 15. 2[) and the analytically computed mean escape time found 
using the fifth-order center manifold equation [Eqs. (|A.4a[) - (|A.4dp ]. Figure l5T47 a) also 
contains lines of best fit passing through the data points associated with numerical 
and analytical computation. 

One can see from Fig. I5.4f a) that there is good agreement between the two meth- 
ods. The slope of the line of best fit through the numerical data is m num = 0.2536, 
while the slope of the line of best fit through the analytical data is m ana = AU = 
0.2591. One also can see in Fig. 15.4( a) that there is a slight discrepancy between 
the analytical and numerical y-intercept values. This is due to the fact that in the 
numerical computation, the escape time was based on the time for the particle to 
cross the barrier and descend partially into the second basin of attraction, while the 
analytical escape time was based on the time for the particle to reach the top of the 
barrier. 

The two methods continue to agree well as e is increased to e = 0.1 [Fig. I5.4f b)] 
and e = 0.14 [Fig. I5.4f c)]. Beyond e = 0.14, the analytical result begins to diverge 
from the numerical result, and this divergence increases as e increases. An example 
of this divergence is shown for e = 0.2 in Fig. I5.4f d). The slopes of the lines of 
best fit are as follows: for e = 0.1, m num = 0.2505 and m ma = AU = 0.2624; for 
e = 0.14, m num = 0.2490 and m ana = AU = 0.2385; and for e = 0.2, m nura = 
0.2466 and m ana = AU — 0.1859. Note that by linearizing the deterministic form of 
Eqs. (|3.1aj) - (|3.1bj) about the stable, attracting point (1,0), one finds that the critical 
damping value is e cr = 0.125. This value agrees well with the value at which the good 
comparison between the analytical result and numerical result begins to break down. 

As stated in Sec. 13.31 we have computed the analytical escape rate of the particle 
using the stochastic center manifold. Again, the stochastic correction is minimal, and 
therefore these analytical mean escape times are very close in value to those found 
using the "naive" approach. We do not show figures comparing this analytical result 
with the numerical result since there is no noticeable difference from the plots shown 
in Figs. IE47a)-ICTd). 

3.5. Numerical Computation of Escape Prehistory. For each of the 10, 000 

particles that were initially placed in one of the attracting basins and which later 
escaped from this basin, across the saddle, and into the other basin of attraction, we 
retain t = 200 worth of the particle's path prior to escape. By creating a histogram 
representing the probability density, p^, of this escape prehistory [14] . one can see 
which regions of the phase space are associated with a high or low probability of 
particle escape. 

An example of a histogram of escape path prehistory is shown in Fig. I5.5f a) for 
e = 0.1 and a = 0.3 (so that D = a 2 /2 = 0.045). The color-bar values of Fig. l53ta) 
have been normalized by 10 5 . The threshold value in Fig. I5.5f a) is actually about 
9000. Therefore, any histogram box containing less than 9000 events shows up as 
white on the histogram. Figure [575T b) is the same as Fig. I5.5f a). but with adjusted 
color-bar values. Doing this enables one to obtain a clearer view of the escape path 
prehistory along the separatrix and near the saddle. Figure 15.5T b) also includes the 
escape path prehistory for one particular particle. To avoid clutter, we have overlaid 
only t = 50 of path prehistory for this one particle. 

Figure I5.6f a) shows Fig. I5.5f a) overlaid with the graph of the slow manifold 
equation, y — x — x 3 . Even though this equation is found by setting e = in 
the deterministic version of Eqs. p. lap and (|3.1b[) . we see in Fig. 15.6( a) that the slow 
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manifold lies very close to the region associated with the highest probability of escape. 

Similarly, Fig. l5.6f b) shows Fig. 15.5( a) overlaid with the graphs of the third-order, 
fourth-order, and fifth-order center manifold equations. Each of these equations may 
be found by including terms of the appropriate order from Eq. ()3.5a|) . As in Fig. l5.5f a). 
the color-bar values of Figs. 15.61 a) and !5.6f b) have been normalized by 10 5 . 

One can see from Fig. I5.6f b) that the third-order and fourth-order center man- 
ifolds essentially bound the entire region of escape path prehistory, while the fifth- 
order manifold lies along the region of highest probability of escape. Although it is 
not shown, it should be noted that the optimal escape path [Eq. (|3.8[) ] associated with 
the third-order center manifold is a heteroclinic orbit from x = \/l — e + 2e 2 to x = 
that lies directly on top of a section of the third-order center manifold [solid, green 
line in Fig. I5.6f b)]. Similarly, the optimal escape paths associated with higher order 
center manifolds lie directly on top of a section of the corresponding center manifold. 

Additionally, one could overlay the histogram of escape path prehistory with the 
average stochastic center manifold given by Eq. (|3 . 24|) . However, since the stochas- 
tic correction appears at order C(e 4 cr 2 ), there is no noticeable difference from the 
manifolds shown in Figs. I5.6f a) and l5.6f b). Therefore, plots of the average stochastic 
manifold are not shown. 

4. Conclusions. A general procedure consisting of elements taken from deter- 
ministic and stochastic manifold theory and large fluctuation theory is developed and 
used to understand the underlying structure of a stochastic dynamical system with 
two well-separated time scales. 

As a first step towards this goal, we have applied the procedure to a generic 
2D, singularly perturbed, damped, Duffing oscillator system with additive, Gaussian 
noise. The deterministic center manifold equation is found by neglecting the stochastic 
terms. By "naively" adding a noise term to the equation that describes the dynamics 
on the center manifold, the path integral formalism of large fluctuation theory can 
be used to analytically compute the optimal escape path of the particle from one 
basin of attraction to another basin of attraction as well as the particle's escape 
rate. Comparison of the analytical result with numerical computations shows excellent 
agreement if the system is in the over-damped regime. There is a shift due to the pre- 
factor of the distribution, but the exponent is quite accurate. While it is possible to 
correctly determine the pre-factor using spectral methods, the technique works only 
for ID problems [38j . When the system enters the under-damped regime, we must 
consider the 2D topology rather than the ID manifold topology that is associated with 
the over-damped regime [T7] . We are currently performing the analysis associated with 
the under-damped regime, and this will be presented elsewhere. 

Additionally, we used numerical computation to create a histogram of the escape 
path prehistory distribution. The histogram enables one to see regions of the phase 
space that are associated with high and low probability of particle escape from the 
basin of attraction. By overlaying the histogram with the deterministic slow manifold, 
we see that this manifold lies close to the region associated with the highest probability 
of escape. By comparing the histogram with deterministic center manifolds of various 
orders, we see that the third-order and fourth-order center manifolds essentially bound 
the entire region of escape, while the fifth-order center manifold lies very close to 
the region associated with the highest probability of escape. Therefore one can use 
these manifolds to accurately describe the location of escape regions. The example 
considered here is that of a quartic potential which is a generic potential for a pitchfork 
bifurcation. Therefore, for systems whose center manifolds yield pitchfork dynamics, 
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we expect similar results to hold when additive noise is considered. 

Knowledge of the location of the regions of high and low probability of escape 
(whether from numerical or analytical results) is extremely useful if one wishes to 
optimize the amount of time spent in a particular region of phase space. For example, 
one might wish to keep an autonomous glider in some basin of attraction. Instead of 
constantly actuating the glider controls to keep the glider stationed at a particular 
location, one can station the glider initially in a region of low probability of escape. 
Then, if the glider enters a region of high probability of escape, one actuates the 
controls to move the glider back to a region with low probability of escape and the 
controls are switched off. This is similar to the stochastic control used in [451 H3] ■ 

As an example, consider a particle whose dynamics are described by Eqs. (|3. la|) 
and (|3.1b[) . Figurc l5~?T a') shows the location of the particle (x coordinate) as a function 
of time t. One can see in Fig. I5.7f a) that the particle moves around the basin of 
attraction until eventually the stochastic effects cause the particle to escape from the 
basin at t w 1290. 

By implementing a simple control that pushes the particle into a region of phase 
space that has a low probability of escape whenever the particle enters a region of 
high probability of escape, it is possible to keep the particle in the basin of attraction 
for a much longer amount of time. Figure I5.7f b) shows the location of the particle 
(x coordinate) as a function of time t with the application of control. One can see 
in Fig. I5.7f b) that control allows one to keep the particle in the basin of attraction 
until t = 10, 000, when the simulation was stopped. Preliminary results show that the 
application of control increases the mean time to escape from the basin and decreases 
the total number of particles which escape in a given period of time. A complete 
analysis of the application of control will appear elsewhere. 

Even though we showed that analytical results using a "naive" approach agree 
very well with the numerical results in the over-damped regime, to more accurately 
describe the stochastic effects, we derived the normal form coordinate transform. 
The normal form enabled us to find the stochastic center manifold. However, the 
stochastic effects appeared in the manifold equation at high order, so this correction 
has a minimal effect. Therefore, for this problem at least, one can use the simpler 
and much less time consuming "naive" approach. 

It should be noted that there are systems where one should not rely on the 
"naive" approach. For example, in a Susceptible-Exposed-Infected-Recovered (SEIR) 
epidemiological model, there are terms at low order in the normal form transform 
which cause a significant difference between the average stochastic center manifold 
and the deterministic manifold j^T]. Therefore, when working with the SEIR model, 
one must use the stochastic normal form coordinate transform approach to obtain the 
correct projection of the noise onto the center manifold. 

Figure 15.81( a) compares the fraction of the population that is infected with a 
disease, /, computed using the complete, stochastic system of equations of the SEIR 
model with the time series of / computed using the reduced system of equations of 
the SEIR model that is based on the deterministic center manifold with a "naive" 
replacement of the noise terms. One can see that the solution computed using the 
reduced system incorrectly predicts the time and amplitude of the initial outbreak 
and quickly becomes out of phase with the solution of the complete system. Although 
not shown, the poor agreement, in phase and amplitude, between the two solutions 
continues for long periods of time. 

On the other hand, Fig. 15.8T b) compares the time series of / computed using 
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the complete, stochastic system of equations of the SEIR model with the time series 
of / computed using the reduced system of equations of the SEIR model that is 
found using the stochastic normal form coordinate transform. There is excellent 
agreement between the two solutions. The initial outbreak is successfully captured 
by the reduced system, and the reduced system continues to accurately predict the 
phase and amplitude of outbreaks over long periods of time. 

The previous general analysis has been performed only for the 2D singularly 
perturbed Duffing system. Of interest is the application of the theory to more realistic 
stochastic dynamical systems. Beyond that, we plan to apply the theory to fully 3D 
systems as well as to actual occanographic data. 
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Appendix A. Escape Rate Using Stochastic ("Naive") Center Mani- 
folds. Using the fourth-order stochastic ("naive") center manifold dynamical equa- 
tion given by 

(A.l) x = x - x 3 - ex + 2e 2 x - 5e 3 x + 4ex 3 + V2D(j)(t), 

one finds that the escape rate from the attractor located at 

x = x a = v /(l-e + 2e 2 -5e 3 )/(l-4e) 
to the saddle located at x = x s = is given as follows: 

(A.2.) ^,). igl!W ^ m 
where 

(A.2b) U"(0) = -1 + e - 2e 2 + 5e 3 , 

(A.2c) U"(x a ) = 2 - 2e + 4e 2 - 10e 3 , 

AU = |( -1 + 6e - 13e 2 + 34e 3 - 70e 4 + 76e 5 
(A.2d) -105e 6 + 100e 7 ) / [4 (l - 8e + 16e 2 )] | . 

Using the fifth-order stochastic ("naive") center manifold dynamical equation 
given by 

(A.3) x = x - x 3 - ex + 2e 2 x - 5e 3 x + 4ex 3 + 14e 4 .T - 20e 2 x 3 + V2Dcj)(t), 
one finds that the escape rate from the attractor located at 

x = x a = y/(l-e + 2e 2 - 5e 3 + 14e 4 )/(l - 4e + 20e 2 ) 
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to the saddle located at x = x s = is given as follows: 

(A.4a) W(e,D) = ^^M exp (-AU/D), 
where 

(A.4b) U"(0) = -l + e-2e 2 + 5e 3 - 14e 4 , 

(A.4c) U"(x a ) = 2 - 2e + 4e 2 - 10e 3 + 42e 4 , 

Af7 = |( 1 - 2e + 5e 2 - 14e 3 + 42e 4 - 48e 5 +81e 6 - 140e 7 + 196e 8 ) 
(A.4d) / [4(-l + 4e-20e 2 )]|. 

Appendix B. Second Iteration Details. For this second iteration, we seek 
a correction to the x coordinate (slow process) with the form 

(B.la) x = X + £(t,X,Y) + ..., 

(B.lb) X' = F(t,X,Y) + ..., 

where £ and F are small corrections. Substitution of Eqs. (|3.17a[) - (jB.lb[) into Eq. 
(|3.12a|) leads to 

<-> * + £ + £% + g»-<r + x-*) + -* 

Replacing X' = dX/dr with F [Eq. |B~Tb|) ]. replacing dY/dr with —Y [Eq. (|3.17bj) ]. 
and ignoring the term d^/dX ■ F since it is a product of small corrections gives the 
following equation: 

(B.3) F + ?f-Y§ 7 = e(Y + X-X 3 ) + ea0. 

or oY 

Equation (|B.3[) must now be solved for F and £. As in the first step, we employ 
principle (UJ) and keep the evolution equation [Eq. (|B.lb[) ] as simple as possible. How- 
ever, since the terms e(X — X 3 ) located on the right-hand side of Eq. (|B.3[) do not 
contain r or Y, these terms must be included in F. Therefore, one piece of F will be 
F = e(X - X 3 ). 

The remaining deterministic term on the right-hand side of Eq. (|B.3[) contains Y. 
This term can therefore be integrated into £. The equation to be solved is 

(B-4) ~ Y W = €Y > 

whose solution is given as £ = — eY. 

To abide by principle we would like to integrate the stochastic piece on the 
right-hand side of Eq. I|B.3[) into £, by solving the equation 

(B.5) dt/dr = ea<j). 
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However, the solution of Eq. (|B.5|) is given by 
(B.6) £ = «j / <f>dr, 



which has secular growth like a Wiener process. Since this would violate principle J!}, 
we must let F = ta<f>. 

Putting the three pieces together yields £ = — eY and F = e(X — X 3 ) + eacj). 
Therefore, the new approximation of the coordinate transform and its dynamics are 
given by 

(B.7a) x = X -eY + 0(( 3 ), 

(B.7b) X' = e(X-X 3 ) + ecrcf) + 0(C 3 )- 

Appendix C. Third Iteration Details. To find the corrections, 77 and G, 
we substitute Eqs. (|3.14a[) and (|3.14b|) into Eq. (|3.12b|) . which leads to the following- 
equation: 

„ drj di] dX dr] , 

Substitution of the specific form of x given by Eq. (|B.7a[) , the specific form of X' = 
dX/dr given by Eq. (|B.7b[) . and substitution of 77 = X — X 3 and G — [see 
Eqs. (j3~T7a| and (|3T7b|) ] into Eq. (fUTT]) gives one the following evolution equation 
driven by the updated residual of Eq. ()3. 12b[) : 

G + ^ - Y^- +7]=e (-X —Y + AX 3 — 3X 5 + 3X 2 Y) 
or oY y ' 

(C.2) + ea (-(j) + iX 2 4>) + e 2 (-3XY 2 ) + e 3 Y 3 . 

We first consider the deterministic terms on the right-hand side of Eq. (|C.2[) . 
Principle ((4]) is employed to keep the evolution equation [Eq. (|3.14b[) ] as simple as 
possible, and since the terms — eX, 4eX 3 , — 3eX 5 arc not functions of r or Y, we let 
7/ = e(— X + AX 3 — 3X 5 ). Consideration of the term — 3e 2 XY 2 leads one to solve the 
following equation: 

(C.3) v - Y ^L = -3e 2 XY 2 . 

The solution of Eq. jC3| is rj = 3e 2 XY 2 . 

The last deterministic terms on the right-hand side of Eq. (|CJ.2[) arc — eY and 
3eX 2 Y. Since these two terms can't be integrated into 77, they are included in G. The 
higher order term e 3 Y 3 will be ignored until a later iteration. 

We now consider the stochastic terms on the right-hand side of Eq. (|C.2[) . The 
equation to be solved is 

(C.4) + ?; = -&r<f> + 'ieaX 2 (i>, 

OT 

whose solution is given by 

r r 

(C.5) 77=-ecr / cxp [— (t - s)]4>{s) ds + ieaX 2 \ exp [-(r - s)](j>(s) ds. 
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If wc define 

T 

(C.6) J exp[-(r-s)]0(s)ds = e" r 

— oo 

then Eq. (|C.5[) can be written as follows: 

(C.7) r) = -eaer T * + 3ecrX 2 e~ T * 4>. 

Putting all of the r\ and G pieces from this third iteration together leads to the 
updated approximation of the coordinate transform and its evolution equation given 
by Eqs. (|338a|) and (f3~18bl) . 

Appendix D. Fourth Iteration Details. To find the corrections, £ and F, 
we substitute Eqs. (jB.lap and (|B.lb[) into Eq. (|3. 12a|) . which leads to the following 
equation: 

(m) F+ 0-r + 0X^ + dY^ = ey + e ^ 

Substitution of the specific forms of y [Eq. (|3.18aj) ]. Y' = dY/dr [Eq. (|3.18b|) ]. £, 
[Eq. (|B.7aj) ]. and F [Eq. (|B.7b|) ] into Eq. (jD.lj) leads to the following evolution equa- 
tion driven by the updated residual of Eq. (|3.12a[) : 

F+^-- Y^X =e 2 (-X + AX 3 - 3X 5 - Y + 3X 2 Y) 
or oY v ' 

(D.2) + e 2 a (-e- T * + 3A 2 e~ T * 4>) + 3e 3 XY 2 . 

It is straightforward (and similar to what has been done in the previous iterations) 
to integrate the deterministic terms on the right-hand side of Eq. (|D.2[) into F and £. 

Consideration of the stochastic terms on the right-hand side of Eq. (|D.2j) means 
that we must solve the following equation: 

(D.3) F+ = e 2 a(-e- T * + 3A 2 e" T * cf>) . 

OT 

As in the second iteration, we can't integrate <j> into £, since this would generate 
secular growth in violation of principle ([1]). Employing principle j5]) to avoid a fast 
time convolution (memory integral) in the slow evolution F leads us to perform an 
integration by parts on each of the terms on the right-hand side of Eq. (|D.3j) so that 

(D.4) F + -^=e 2 a [-4 + e_T * <t>' + (<t> - e ~ T * </»')] • 

It is clear that we now have F = e 2 cr(— </> + 3A 2 0) and £ = e 2 a(e~ T * </> — 3X 2 e~ T * 4>). 

Putting all of the £ and F pieces from this fourth iteration together leads to the 
updated approximation of the coordinate transform and its evolution equation given 
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by Eqs. (jXTOaj) and (|3.19bj) . 

Appendix E. Normal Form Coordinate Transform. 

y =Y + X - X 3 + e (-X + AX 3 - 3X 5 ) + ea (-e~ T * + 3X 2 e~ T * 0) 
+ e 2 (2X - 20X 3 + 42X 5 + 3XY 2 ) 

+ e 2 a (2e- T * <j) + e~ T * e~ r * <j> - 18X 2 e~ T * - 12X 2 e~ T * e~ r * </> 

+ 24X 4 e~ r * (j) +15X 4 e' T * e" r * 0) 
+ e 3 (-Jf + 16X 3 - 66X 5 + 96X 7 - 45X 9 -9XF 2 - Y 3 /2 + 33X 3 Y 2 ) 
+ e 3 a (-e- T * 4> + 18X 3 re" T * - 6Xy e ~ r * + 15X 2 e~ r * <j> 

+ 45X 6 e _T * - 51X V T * (/> + 18X 6 e~ T * e~ T * (j) 

- 2AX i e~ T * e" r * + 6X 2 e~ T * e _T * +3Y 2 e +T * 0) 
+ e 4 (18X 5 Y 2 - 6X 3 Y 2 + 3Y 3 /2 -9Y 3 X 2 /2) 

+ e 4 cr (6XYe~ T * <j> + 54X 5 Ye~ r * - 36X 3 r e " T * + 9X 2 r 2 e +r * 

- 3Y 2 e +T * + 3r 2 e +r * e~ T * 4> -9X 2 Y 2 e +T * e~ T * 4>) 

+ eV (-3Xe~ r * (e~ T * 0) 2 - 27X 5 e~ T * (e~ T * 0) 2 
(E.la) +18X 3 e- T * (e- T *0) 2 ) +C(C 4 ), 




x =X - eY + e 2 (Y - 3X 2 Y) + e 2 a (e~ T * - 3X 2 e~ T * <f) 
+ e 3 (6X 2 Y - 12X A Y -2Y - 3XY 2 /2) 

+ e 4 (~9X 6 Y + 3X 4 Y - 3X 2 Y + Y + Y 3 /6 +9XY 2 /2 - 33X 3 Y 2 /2) 
+ e 5 (-Y 3 /2 + 3X 3 Y 2 + 3X 2 Y 3 /2 -9X 5 Y 2 ) 

+ t 3 a (-3e~ T *4>- 33X 4 e~ T * + 24X 2 e~ T * <j) - 15X 4 e~ T * e~ T * <p 

- 6XYe +T * cf> - e~ T * e~ T * <f> -12X 2 e" r * e~ T * 4>) 
+ e 4 cr (e~ T * 4> - 15X 2 e~ r * - 6X 2 e~ T * e~ T * + 51X 4 e" T * 



+ e 5 cr {~MX 3 Ye +T * 4> + 9XYe +T * + 81x 5 Ye +T * <f> - 3X1^ * 
- 27X 5 Ye~ T *<t>+ l%X 3 Ye~ T * - 9X 2 Y 2 e +2T * e +r * </> 
+ 3 (r 2 e +2T * e +T * (f) 12 -3 (Y 2 e +2r * e~ T * 0) /2) 

+ e 6 er (-6Xy e +T * + 54X 3 Fe +T * - 162X 5 Fe +T * +162X 7 Fe +T * 0) 




(E.lb) 



- Y + e (-Y + 3X 2 Y) + e 2 (Y - 6X 2 Y + 9X 4 y) 
+ e 3 a (6XY(f) - 18X 3 Ycj)) 

+ e 4 a (-6XYcj) - hAX 5 Y<t> +36X 3 Y<t>) + (D{( 4 ), 



- 45X 6 e~ T *4> + 24X 4 e~ T * e~ T * cj> - 18X 6 e" r * e~ T * 4> 

+ 3XYe +T * + 3XYe~ T * - 9X 3 Ye +T * - 9X 3 Ye~ T * (f> 

-3Y 2 e +2T * e +T * 0) 
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+ eV [3X (e- T *tf>) /2- 9X 3 (e~ T * 4>) + 27 X 5 (e^ 7 * </>) /2 
+ 3Xe~ T * (e~ T * cj)) 2 - 18X 3 e~ T * (e~ r * (ff 
(E.lc) +27X 5 e~ T * ( e - r *</,) 2 ) +C(C 4 ), 



X' =e (X - X 3 ) + e<7(f> + e 2 (-X + iX 3 - 3X 5 ) + e 2 a (-0 + iX 2 qb) 
+ e 3 (2X - 20X 3 + 42X 5 ) 
+ e 4 {-X + l&X 3 - 66X 5 96X 7 - 45X 9 ) 

+ e 3 cr (30 - 2AX 2 4> + 33X 4 0) + eV (-0 + 15X 2 - 51X 4 c/> + 45X 6 9 
+ eV (~6X(j)er T * + 18X 3 0e~ T * (f>) 
(E.ld) + eV 2 (-9X0e" r * <p + bAX 3 (jye~ T * -81X 5 0e" r * <p) + C(C 4 ), 



where 

(E.2) e +T *</>= / exp [(r - s)]<f>(s) ds. 



REFERENCES 

[1] L. Arnold, Random Dynamical Systems, Springer- Verlag, 1998. 

[2] L. Arnold and P. Imkeller, Normal forms for stochastic differential equations, Probab. 
Theory Rel., 110 (1998), pp. 559-588. 

[3] N. Berglund and B. Gentz, Geometric singular perturbation theory for stochastic differential 
equations, J. Differ. Equations, 191 (2003), pp. 1-54. 

[4] L. Billings, E. M. Bollt, and I. B. Schwartz, Phase-space transport of stochastic chaos in 
population dynamics of virus spread, Phys. Rev. Lett., 88 (2002), paper 234101. 

[5] E. M. Bollt, L. Billings, and I. B. Schwartz, A manifold independent approach to under- 
standing transport in stochastic dynamical systems, Physica D, 173 (2002), pp. 153—177. 

[6] P. Boxler, A stochastic version of center manifold theory, Probab. Theory Rel., 83 (1989), 
pp. 509-545. 

[7] J. Carr, Applications of Centre Manifold Theory, Springer- Verlag, 1981. 

[8] T. Carr, L. Billings, I. Schwartz, and I. Triandaf, Bi-instability and the global rate of 

unstable resonant orbits in a driven laser, Physica D, 147 (2000), pp. 59-82. 
[9] P. H. Coullet, C. Elphick, and E. TlRAPEGUl, Normal form of a H op f bifurcation with noise, 
Phys. Lett. A, 111 (1985), pp. 277-282. 
[10] R. E. L. DeVille, E. Vanden-Eijnden, and C. B. Muratov, Two distinct mechanisms of 
coherence in randomly perturbed dynamical systems, Phys. Rev. E, 72 (2005), paper 031105. 
[11] M. I. Dykman. Private communication. 

[12] M. I. Dykman, Large fluctuations and fiuctuational transitions in systems driven by coloured 

gaussian noise: A high-frequency noise, Phys. Rev. A, 42 (1990), pp. 2020-2029. 
[13] M. I. Dykman and B. Golding, Controlling large fluctuations: Theory and experiment, in 

Stochastic Processes in Physics, Chemistry, and Biology, J. A. Frcund and T. Poschel, eds., 

Springer- Verlag, 2000, pp. 365-377. 
[14] M. I. Dykman, P. V. E. McClintock, V. N. Smelyanski, N. D. Stein, and N. G. Stocks, 

Optimal paths and the prehistory problem for large fluctuations in noise-driven systems, 

Phys. Rev. Lett., 68 (1992), pp. 2718-2721. 
[15] M. I. Dykman, E. Mori, J. ROSS, and P. M. Hunt, Large fluctuations and optimal paths in 

chemical-kinetics, J. Chem. Phys., 100 (1994), pp. 5735-5750. 
[16] M. I. Dykman, I. B. Schwartz, and A. S. Landsman, Disease extinction in the presence of 

random vaccination, Phys. Rev. Lett., 101 (2008), paper 078101. 



ESCAPE RATES IN A STOCHASTIC ENVIRONMENT 



23 



[17] M. I. Dykman, I. B. Schwartz, and M. Shapiro, Scaling in activated escape of underdamped 
systems, Phys. Rev. E, 72 (2005), paper 021102. 

[18] S. J. B. ElNCHCOMB AND A. J. McKane, Using path-integral methods to calculate noise-induced 
escape rates in bistable systems: The case of quasi-monochromatic noise, in Fluctuations 
and Order: The New Synthesis, Mark Millonas, ed., Springer- Verlag, 1996, pp. 139—154. 

[19] C. C. Eriksen, T. J. Osse, R. D. Light, T. Wen, T. W. Lehman, P. L. Sabin, J. W. 

Ballard, and A. M Chiodi, Seaglider: A long-range autonomous underwater vehicle for 
oceanographic research, IEEE J. Oceanic Eng., 26 (2001), pp. 424-436. 

[20] R. P. Feynman and A. R. HlBBS, Quantum Mechanics and Path Integrals, McGraw-Hill, Inc., 
1965. 

[21] E. Forgoston, L. Billings, and I. B. Schwartz, Accurate time series prediction in 
reduced stochastic epidemic models, 2009. submitted, available at preprint archive: 
[http: / /arxiv.org/abs/0903. 1038 1 

[22] M. I. Freidlin, On stable oscillations and equilibriums induced by small noise, J. Stat. Phys., 
103 (2001), pp. 283-300. 

[23] M. I. Freidlin and A. D. Wentzell, Random Perturbations of Dynamical Systems, Springer- 
Vcrlag, 1984. 

[24] L. Gammaitoni, P. Hanggi, P. Jung, and F. Marchesoni, Stochastic resonance, Rev. Mod. 
Phys., 70 (1998), pp. 223-287. 

[25] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural 
Sciences, Springer- Verlag, 2004. 

[26] R. Graham and T. Tel, Existence of a potential for dissipative dynamical systems, Phys. Rev. 
Lett., 52 (1984), pp. 9-12. 

[27] J. Hales, A. Zhukov, R. Roy, and M. I. Dykman, Dynamics of activated escape and its 
observation in a semiconductor laser, Phys. Rev. Lett., 85 (2000), pp. 78—81. 

[28] A. Hamm, T. Tel, and R. Graham, Noise-induced attractor explosions near tangent bifurca- 
tions, Physics Letters A, 185 (1994), pp. 313-320. 

[29] J. A. Hansen and C. Penland, Efficient approximate techniques for integrating stochastic 
differential equations, Mon. Weather Rev., 134 (2006), pp. 3006-3014. 

[30] B. Hauschildt, N. B. Janson, A. Balanov, and E. Scholl, Noise-induced cooperative dy- 
namics and its control in coupled neuron models, Phys. Rev. E, 74 (2006), paper 051906. 

[31] G. Hu, Stationary solution of master- equations in the large-system-size limit, Phys. Rev. A, 
36 (1987), pp. 5782-5790. 

[32] S. K. Hwang, J. B. Gao, and J. M. Liu, Noise-induced chaos in an optically injected semi- 
conductor laser model, Phys. Rev. E, 61 (2000), pp. 5162-5170. 

[33] A. Kamenev and B. MeersON, Extinction of an infectious disease: A large fluctuation in a 
nonequilibrium system, Phys. Rev. E, 77 (2008), paper 061107. 

[34] E. Knobloch and K. A. Wiesenfeld, Bifurcations in fluctuating systems: The center- 
manifold approach, J. Stat. Phys., 33 (1983), pp. 611-637. 

[35] Y-C Lai, Z. Liu, L. Billings, and I. B. Schwartz, Noise-induced unstable dimension vari- 
ability and transition to chaos in random dynamical systems, Phys. Rev. E, 67 (2003), 
paper 026210. 

[36] D. G. Luchinsky, P. V. E. McClintock, and M. I. Dykman, Analogue studies of nonlinear 

systems, Rep. Prog. Phys., 61 (1998), pp. 889-997. 
[37] R. S. Maier and D. L. Stein, Escape problem for irreversible systems, Phys. Rev. E, 48 (1993), 

pp. 931-938. 
[38] B. Meerson. Private communication. 

[39] M. Millonas, ed., Fluctuations and Order: The New Synthesis, Springer- Verlag, 1996. 

[40] N. S. NAMACHCHIVAYA, Stochastic bifurcation, Appl. Math. Comput., 38 (1990), pp. 101-159. 

[41] N. S. NAMACHCHIVAYA AND Y. K. Lin, Method of stochastic normal forms, Int. J. Nonlinear 

Mech., 26 (1991), pp. 931-943. 
[42] A. J. Roberts, Normal form transforms separate slow and fast modes in stochastic dynamical 

systems, Physica A, 387 (2008), pp. 12-38. 
[43] W. RijMELlN, Numerical treatment of stochastic differential equations, SIAM J. Numer. Anal., 

19 (1982), pp. 604-613. 

[44] Z. SCHUSS AND A. Spivak, Where is the exit point?, Chem. Phys., 235 (1998), pp. 227-242. 
[45] I. B. Schwartz, L. Billings, and E. M. Bollt, Dynamical epidemic suppression using 

stochastic prediction and control, Phys. Rev. E, 70 (2004), paper 046220. 
[46] I. B. Schwartz, L. Billings, M. Dykman, and A. Landsman, Predicting extinction rates in 

stochastic epidemic models, J. Stat. Mech. -Theory E., (2009), paper P01005. 
[47] J. Sherman, R. E. Davis, W. B. Owens, and J. Valdes, The autonomous underwater glider 

'Spray', IEEE J. Oceanic Eng., 26 (2001), pp. 437-446. 



E. FORGOSTON AND IRA B. SCHWARTZ 



D. C. Webb, P. J. Simonettt, and C. P. Jones, Slocum: An underwater glider propelled by 
environmental energy, IEEE J. Oceanic Eng., 26 (2001), pp. 447-452. 

A.D. Wentzell, Rough limit theorems on large deviations for Markov stochastic processes, I, 
Theor. Probab. Appl., 21 (1976), pp. 227-242. 




Fig. 5.1. (a) (j> op t(t) with D = 0.05 for e = 0.02 (red), e = 0.1 (blue), e = 0.25 (green), e = 0.5 
(cyan), t = 1.0 (magenta), and e = 1.5 (black), (b) A close-up view of a section of Fig. [5. lY a). 




Fig. 5.2. Contours of numerically computed mean escape times (plotted as the natural log of 
the mean escape time) as a function of e and 1/D for 10,000 particles. 
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Fig. 5.3. (a) Vertical slice of Fig. \5.S\ taken at e = 0.1. (b) Vertical slices of Fig. \5.2\ taken 
at e = 0.1 (black, "circle" markers), e = 0.2 (red, "square" markers), e = 0.3 (blue, "triangle" 
markers), e = 0.4 (green, "cross" markers), e = 0.5 (magenta, "diamond" markers), and e = 1.0 
(cyan, "asterisk" markers). The data points in both (a) and (b) are overlaid by a line of best fit. 
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Fig. 5.4. Comparison of numerical data ("circle" markers) and linear fit (solid line) with 
analytical data ("square" markers) and linear fit (dashed line) of the natural log of the mean escape 
time as a function of 1/D for (a) e = 0.02, (b) e = 0.1, (c) t = 0.14, and (d) e = 0.2. 
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Fig. 5.5. (a) Histogram of escape path prehistory (fort = 200 of prehistory) of 10, 000 particles 
with e = 0.1 and a = 0.3. The color-bar values have been normalized by 10 s , and the threshold 
value is about 9000. (b) Same as Fig. 15. 3( a), but with adjusted color-bar values, and including one 
particular particle's escape path prehistory (showing only t = 50 of prehistory). 
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Fig. 5.6. Escape path prehistory histogram of Fia. Wlft a) overlaid with (a) the graph of the slow 
manifold eguation, and (b) the graphs of the third-order (solid, green line), fourth-order (dashed, 
green line), and fifth- order (solid, red line) center manifold eguations given by Eq. i3. 5a\) . For both 
Figs. U)7bY a) and \5.bt b), the color-bar values have been normalized by 10 5 , and the threshold value 
is about 9000. 
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Fig. 5.7. Location of a particle (x coordinate) as a function of timet (a) without control, 
(b) with control. 
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Fig. 5.8. Time series of the fraction of the population that is infected with a disease, I, 
computed using the complete, stochastic system of equations of the SEIR model (red, solid line), 
and (a) computed using the reduced system of equations of the SEIR model that is based on the 
deterministic center manifold with a "name" replacement of the noise terms (blue, dashed line), 
and (b) computed using the reduced system of equations of the SEIR model that is found using the 
stochastic normal form coordinate transform (blue, dashed line). 



